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The distribution of entanglement between 
macroscopically separated parties represents a 
crucial protocol for future quantum information 
networks pQ[2]. Surprisingly, it has been theoret- 
ically shown that two distant systems can be en- 
tangled by sending a third mediating system that 
is not entangled with either of them [S[ [4] . Such a 
possibility seems to contradict the intuition that 
to distribute entanglement, the transmitted sys- 
tem always needs to be entangled with the sender. 
Here, we experimentally distribute entanglement 
by exchanging a subsystem and successfully prove 
that this subsystem is not entangled with either 
of the two parties. Our implementation relies on 
the preparation of a specific three-mode Gaus- 
sian state containing thermal noise that demol- 
ishes the entanglement in two of the three bipar- 
tite splittings. After transmission of a separa- 
ble mode this noise can be removed by quantum 
interference. Our work demonstrates an unex- 
pected variant of entanglement distribution and 
improves the understanding necessary to engineer 
multipartite quantum information networks. 

The principle of entanglement distribution by separa- 
ble states is illustrated in Fig. [I] In the beginning of the 
protocol Alice possesses two separable modes A and C, 
while Bob possesses the mode B, which is also separa- 
ble from Alice's modes. In a first step Alice sends the 
ancilla mode C, which is neither entangled with mode A 
nor with mode B, to Bob. To obtain two-mode entangle- 
ment Bob mixes his modes B and C in the second step of 
the protocol. One output mode is then discarded, while 
the other one turns out to be entangled with A. The dis- 
tributed entanglement can be used for further quantum 
information protocols [5] , such as quantum teleportation 
[7] and quantum key distribution [8] . 

This remarkable and seemingly paradoxical protocol is 
made possible by a specific structure of quantum corre- 
lations in the underlying three- mode Gaussian state. For 
the protocol to work the state must be separable with 
respect to the B|AC and C|AB splittings and insepara- 
ble with respect to the A|BC splitting. According to the 
classification introduced in [9 , we therefore need a three- 
mode Gaussian entangled state belonging to Class III. 

Our protocol thus explores the rich structure of mul- 
timode entanglement, which can exhibit more complex 
properties and features than two-mode entanglement 
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Figure 1: Principle of entanglement distribution by 
separable states. In the beginning Alice possesses the two 
separable modes A and C. Both modes are also separable with 
respect to Bob's mode B. Alice sends mode C to Bob and he 
combines his mode B with the received mode C. Finally, Alice 
and Bob share an entangled system A|B', which can be traced 
back to the initial entanglement for the A|BC splitting. 



and which represents a valuable resource for lots of 
applications ranging from local realism tests [10] to 
one-way quantum computing [TTHT3] . 



Three-mode state preparation 

Our experimental setup is depicted in Fig. [2j The 
initial three-mode Gaussian state is prepared by an 
independent source and is distributed between Alice 
and Bob. The preparation starts with a squeezed state, 
which interferes with a vacuum state at a balanced beam 
splitter. The beam splitter output A is sent directly 
to Alice, while the other output is superimposed with 
a thermal state at a second balanced beam splitter. 
After the state preparation Alice possesses modes A 
and C, while Bob holds mode B. The separability 
properties of this three-mode state (ABC) were checked 
by a tomographic reconstruction of the full three-mode 
covariance matrix with the homodyne detectors HDa, 
HDb and HDc and found to be separable with respect 
to the B|AC and C|BA splittings. In the next step 
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Figure 2: Experimental setup of the entanglement dis- 
tribution by separable states. The three-mode state is 
prepared by overlapping a squeezed state, a vacuum state 
and a thermal state at two balanced beam splitters. After 
modes A and C had been sent to Alice and mode B had been 
sent to Bob, the separability properties were checked with bal- 
anced homodyne detectors HDa, HDb and HDc- Afterwards 
Alice sent the separable mode C to Bob. By overlapping the 
modes C and B at another balanced beam splitter BS3, entan- 
glement was established between Alice and Bob, which was 
verified with homodyne detectors HDa and HD B / . 

Alice sends the mode C to Bob, where the modes B 
and C interfere at a balanced beam splitter with the 
appropriate phase to get rid of the correlated noise. This 
step creates two-mode entanglement between Alice and 
Bob, which is verified by measuring the Duan criterion 

El 

Var(X A - X w ) + Var(P A + Pb>) < 4 (1) 

using the homodyne detectors HDa and HDb'. Here, 
X and P describe the amplitude and phase quadrature 
operators, respectively. They are normalized to the shot 
noise level, i. e. Var(X) = Var(P) = 1 for a vacuum 
state. 

Requirements for the three-mode state 

For investigating the separability properties of the 
three-mode state (ABC) we apply the positive partial 
transposition criterion (PPT) [15j [16] to the measured 
state. This criterion is both necessary and sufficient for 
bipartite splittings of Gaussian states with N modes 
with only a single mode on one side (1\N — 1) [17 . The 
three-mode state is separable with respect to mode k 
if the corresponding covariance matrix of the partially 
transposed state fulfills the uncertainty relation 

7 T(fc) -ffi>0, (2) 
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Figure 3: Separability analysis. Theoretical simulations of 
the PPTc value with respect to the added thermal noise are 
depicted. The solid lines correspond to 10 dB initial squeezing 
with different losses. To obtain the required PPT value (> 1) 
a minimal loss of 33% is necessary. The point of intersection 
with the threshold of 1 is exclusively depending on the loss 
and not on the actual initial squeezing value. The parameters, 
which were used in our measurements, are marked with the 
magenta cross. 

with Q = 0J* =1 J, where J = ^ ^j. This criterion 

is equivalent to finding the symplectic eigenvalues of the 
covariance matrix of the partially transposed state. If 
the smallest symplectic eigenvalue in the following 
called PPT value, is below 1, the state is inseparable with 
respect to the k\ij splitting. For details see supplemental 
material. 

We named the PPT values for the different split- 
tings after the single mode: PPTa for the (A|BC) split- 
ting, PPT B for (B|AC), and PPT C for (C|AB). Our 
protocol thus requires PPTa < 1 (= inseparable) and 
PPTb, PPTc > 1 (= separable) to verify the appropriate 
three-mode state for distributing entanglement by sepa- 
rable states. 

Within the experimental setup we can vary the follow- 
ing critical parameters: the size of the thermal state as 
well as the variances of the squeezed and anti-squeezed 
quadratures of the squeezed state. The latter two can be 
changed independently of each other by variation of the 
pump power of the squeezed light source and by variation 
of additional losses. 

We performed theoretical simulations to analyse the 
influence of the squeezed and thermal states on the sep- 
arability properties of the generated three-mode state. 
Figure [3] shows the PPTc value of the three- mode state 
versus the noise power of the thermal state. Due to the 
symmetry of the setup PPTb and PPTc are identical. 
The PPTa value is not depicted, since this value is al- 
ways below 1 when a squeezed state is used as an input 
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State Aj n . 

The magenta lines show the independence of the ini- 
tial squeezing value on the intersection point with the 
separability threshold unity. Hence, the amount of clas- 
sical noise necessary to demolish the entanglement only 
depends on the optical loss applied to the squeezed state. 
Nevertheless, with higher squeezing values the three- 
mode state is farther pushed into the separable regime 
with respect to the B|AC and C|BA splittings for a suf- 
ficiently large thermal state. Apart from that, higher 
squeezing values also distribute more entanglement to the 
two distant parties. 

From Fig. [3] it is also visible that the entanglement can 
only be demolished by classical noise if the optical loss 
applied to the squeezed state is larger than 33.3%. This 
is exactly the border for which the bipartite entangled 
state, generated by the superposition of the squeezed 
and the vacuum state, is no longer EPR-entangled 
[18]. EPR-entangled states are a subclass of general 
entanglement, exhibiting stronger quantum correlations. 
Indeed the properties of our three-mode state show that 
these correlations are so strong that the entanglement in 
the bipartite splittings cannot be demolished by classical 
noise. 

Experimental Results 

The 21 independent elements of the symmetric 6x6 
three-mode covariance matrix are determined from 
homodyne measurements on modes A, B, and C. For 
each quadrature measurement we recorded 10 6 data 
points. As input states we used a squeezed state with 
-1.8 dB and 5.1 dB noise reduction/amplification in 
the amplitude and phase quadrature, respectively, and 
an elliptical thermal state (hot squeezed state) with 
9.6 dB and 10.2 dB noise amplification. The resulting 
three-mode covariance matrix 7 was measured as 
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This covariance matrix directly leads to the PPT val- 
ues: PPT A = 0.89, PPT B = 1.1 and PPT C = 1.07. 
Thus, the measured state fulfilled the requirements for 
distributing entanglement via separable states. 

Three main effects could in principle cause masking 
the actual presence of entanglement. Two of them can 
also lead to a non- Gaussian state and can thus prohibit 
the application of the separability criterion for Gaussian 
states: phase fluctuations due to imperfect phase locking 
between signal beams and local oscillator beams and the 
generation of the thermal state by random displacements 
of originally squeezed states, where the distribution of 
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Figure 4: Measured PPT values with subtraction of 
detection losses. The magenta curves show the inferred 
PPTb and PPTc values of the measured covariance matrix 7 
for a spectrum of computationally eliminated detection losses. 
Based on independent measurements we estimate the actual 
detection loss to be greater 7% and smaller 22%. These losses 
do not push the PPT values below unity. The successful 
demonstration of our protocol is thus independent of the ques- 
tion whether detection loss should be corrected for or not. In 
contrary, the grey dotted curves, which shows the PPT val- 
ues for another measurement, fall below the threshold, when 
losses > 15% are subtracted and thus do not fulfill the three- 
mode state requirements. 



random displacements can be non-Gaussian. These ef- 
fects are considered in detail in the appendix with the 
result that none of them has any non-negligible effect in 
the presented measurement. 

The third effect is the influence of detection losses. 
Since we are in fact interested in the separability prop- 
erties of the state before homodyne detection, the opti- 
cal loss introduced by the measurement devices has to 
be computationally eliminated. Indeed, the separability 
properties of the state can be altered by a non-perfect 
detection process as depicted in Fig. [4j The blue curves 
represent the PPTb and PPTc values of the covariance 
matrix 7, if optical loss within the homodyne detection is 
subtracted. The black vertical lines mark the regime of 
our estimated detection efficiency (quantum efficiency + 
visibility). We estimate the quantum efficiency of the ho- 
modyne's photodiodes to be about 90%. The visibilities 
of the homodyne detectors were measured before each 
measurement and laid in a regime of 93-98%. For the co- 
variance matrix 7 the detection losses thus are 11% ±5% 
for the homodyne detector HDa, 17% ±5% for HDb and 
16.6% ± 5% for HDc, which leads to an upper bound of 
6% and a lower bound of 22% loss. After subtracting the 
losses from the three-mode covariance matrix the lower 
and upper bounds for the PPT values are 0.85 and 0.87 
for the A|BC splitting, 1.07 and 1.09 for the B|AC split- 
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ting and 1.04 and 1.06 for the C|AB splitting and thus 
fulfill the criteria. This shows the correctness of the sep- 
arability properties regardless whether the detection loss 
is considered to be part of the detected state or not. 

The grey curves in Fig. [4] however represent the PPT 
values of an inappropriate example, where the state is 
inseparable in all splittings, if the estimated detection 
losses are subtracted. Thus, this state does not allow a 
demonstration of entanglement distribution by separable 
states. 

Monte Carlo simulations show that for the 10 6 mea- 
surements per homodyne setting the statistical error bars 
on the symplectic eigenvalues are of the order of 0.001. 
That means, that the inferred separability properties are 
statistically reliable even for the extreme limit of 22% 
loss. 

After the prepared three-mode state had been checked 
for its separability properties, the ancilla mode C, which 
was separable to the modes A and C, was sent to Bob. 
Two-mode entanglement between Alice and Bob was gen- 
erated by superimposing modes C and B at the balanced 
beam splitter BS3 with the appropriate phase, which was 
controlled manually. The Duan criterion resulted in 3.4 
(< 4), which proved that entanglement was successfully 
distributed by separable states. 

In conclusion, we experimentally realized entangle- 
ment distribution by separable states. We showed that 
for this protocol a specific three-mode state is suitable, 
whose thermal noise destroys the entanglement in two 
bipartite splittings, which can later be restored via 
quantum interference. Thereby we could show, that the 
protocol does not work with EPR-entangled states, since 
with states in this class of entanglement separability 
cannot be generated by thermal noise. While the entan- 
glement distribution by separable states seems highly 
counterintuitive in the first place, our protocol provides 
an insight into the underlying physical mechanism be- 
hind this protocol. From a broader perspective our work 
helps to understand the possibilities and restrictions 
offered by mulitmode entangled quantum states and 
future multipartite quantum communication networks. 

Methods 

1. Squeezing Source 

The squeezed states which we used in our protocol were 
generated by parametric amplification. The squeez- 
ing source was a hemilithic optical parametric amplifier 
(OPA) consisting of nonlinear 7% magnesium ox- 
ide doped lithium niobate (MgO:LiNb03) crystal inside 
a standing wave cavity. The cavity was composed of the 
crystal's back surface with a high-reflective coating for 
both wavelengths and the front mirror with a reflectivity 
of 94% for the squeezing field at 1064 nm and a reflectiv- 
ity of 25% for the pump field at 532 nm. These values cor- 
respond to a finesse of 100 for 1064 nm and 4.3 for 532 nm. 
The free spectral range was about 4 GHz due to the cav- 



ity's length of 6.5 mm. The cavity was stabilized by the 
Pound-Drever-Hall technique with a sideband frequency 
of 30 MHz. Therefore, a dim control field of about 2 mW 
entered the cavity from the back end and the length was 
controlled by a piezo electric transducer driving the end 
mirror of the cavity. The crystal's temperature was ac- 
tively stabilized at about 60° C for the phase matching of 
the fundamental and the harmonic field. 

2. Generation of the thermal state 
The thermal state, which we used as the input state Ci n , 
was generated as a classical squeezed state, also called 
a hot squeezed state. For this state preparation we also 
used an optical parametric amplifier (OPA) as described 
above, which reduced the noise in the phase quadrature 
and amplifies noise in the amplitude quadrature. An 
electro-optical modulator applied a random noise on the 
phase quadrature of the vacuum mode entering the OPA. 
The random displacement resulted from the output of a 
homodyne detector measuring shot noise. Therefore, we 
obtained a state with a noise distribution far above the 
vacuum noise in all quadratures. For details see [19] . 
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APPENDIX 

Covariance matrix 

Let Xj and Pj denote the amplitude and phase quadra- 
tures of mode j that satisfy the canonical commuta- 
tion relations [Xj,Pk] = 2iSjk- Here Sjk denotes the 
Kronecker delta and the quadrature operators are nor- 
malized such that the variance of vacuum state quadra- 
tures is equal to 1. When dealing with N modes, it 
is convenient to collect the quadratures into a vector 
Q = (Xi , Pi , • • • , Xj , Pj , • • • , Xn , Pn)- The elements 
of the A^-mode covariance matrix are defined as follows 

ra, 

j jk = ^(AQjAQk + AQ k AQj), (3) 

where AQj = Qj — (Qj) and () denotes averaging over 
a quantum state. The covariance matrix of a physical 
state is symmetric and positive definite and satisfies the 
generalized Heisenberg uncertainty relation [5] 

7 - iJ N > 0, (4) 
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where 

3 = 1 V 7 

is the A/'-mode symplectic form. Symplectic eigenvalues 
of 7 are defined as the positive roots of polynomial 



fluctuations due to imperfect phase locking between sig- 
nal beams and local oscillator beams, and generation 
of hot squeezing by random displacements of originally 
squeezed quadrature where the distribution of random 
displacements may slightly deviate from Gaussian distri- 
bution. In what follows, we will consider each of these 
effects in detail. 



|7 - ihJn\ 



(5) 



where \A\ denotes the determinant of matrix A. The 
symplectic eigenvalues of a covariance matrix of physical 
state satisfy the inequality \i > 1. 

The entanglement and separability of Gaussian states 
may be conveniently analyzed with the help of the partial 
transposition criterion [16J [17] . At the level of quadra- 
ture operators, the partial transposition with respect to 
mode j corresponds to the change of sign of the phase 
quadrature, Pk — » —Pk- The covariance matrix of a state 
partially transposed with respect to mode j thus reads, 



7 



T(k) 



(6) 



where is a diagonal matrix with all diagonal elements 
equal to 1 except for T^k^k — ~ 1- Consider a bipartite 
splitting of an TV-mode Gaussian state with covariance 
matrix 7 such that one party holds mode k and the other 
party possesses all the remaining N—l modes. As shown 
by Werner and Wolf [17] , the partial transposition crite- 
rion is both necessary and sufficient for this case. The 
Gaussian state is entangled with respect to the consid- 
ered bipartite splitting if and only if the lowest symplectic 
eigenvalue nu of satisfies < 1 and otherwise it is 
separable. 



Verification of separability 

In our experiment, we prepare a specific three-mode 
state and we want to demonstrate its separability with 
respect to the bipartite splittings B\AC and C\AB. Since 
the state is generated by mixing single-mode squeezed 
states on beam splitters, it seems natural to assume that 
the state is Gaussian and apply the necessary and suffi- 
cient separability criteria for Gaussian states. However, 
the entanglement of the state may be potentially unde- 
tectable in the Gaussian setting if certain experimental 
imperfections make the state slightly non-Gaussian. 

Here, we identify the main effects that may potentially 
influence the assessment of the entanglement properties 
of the considered state and we demonstrate that the sep- 
arability properties of the experimentally generated state 
remain unchanged when these effects are taken into ac- 
count and compensated for. As illustrated in Fig. S.l, 
the three main effects that may mask the entanglement of 
the state consist of inefficient homodyne detection, phase 



Inefficient homodyne detection 

In the experiment, the three-mode state is probed with 
three balanced homodyne detectors and the covariance 
matrix 7 of the state is determined from the homodyne 
data. The homodyne detectors are characterized by their 
efficiency rjj that is limited by the quantum efficiency of 
the photodiodes, the mode-matching of signal and lo- 
cal oscillator beams and the losses of the passive optical 
elements in the paths of the signal beams. Inefficient 
homodyne detection is equivalent to a combination of a 
transmission through a purely lossy channel with trans- 
mittance 77 followed by ideal homodyne detection. The 
covariance matrix 7^ before transmission through lossy 
channels is then related to the reconstructed covariance 
matrix according to 



7m = S^lS^ + Grj, 



(7) 
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Figure 5: Equivalent scheme of the experimental setup in- 
cluding effects that can potentially mask entanglement of the 
resulting three-mode state. The state preparation begins by 
splitting a single-mode squeezed vacuum on a balanced beam 
splitter. One mode is then mixed with a noisy quasi-thermal 
state produced from a squeezed vacuum state by random dis- 
placements of the initially squeezed quadrature. After mixing 
on beam splitters, each mode can be subject to phase fluctua- 
tions and losses. The purely lossy channel with transmittance 
rjj followed by perfect homodyne detector models a realistic 
homodyne detector with efficiency rjj and the phase fluctu- 
ations account for imperfect phase stability between signal 
and local-oscillator beams in the homodyne detector. Note 
that the phase shift operation commutes with a lossy chan- 
nel hence we can consider the ordering depicted in the figure 
without any loss of generality. 
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and Grj = / — S^S^ where / denotes the identity matrix. 
The covariance matrix 7^ can be determined by inverting 
the linear relation ([7]), 

lL = S- 1 { lM -G n )ST~ 1 . (Q) 

It may happen that the state with covariance matrix jl 
is actually entangled with respect to one or both of the 
splittings B\AC, C\AB although losses induce separa- 
bility of 7m with respect to both these splittings. To 
convincingly verify that this is not the case in our exper- 
iment, we have to take conservative estimates of rjj which 
means that we should consider the lowest possible values 
of rjj consistent with our characterization of homodyne 
detectors and the statistical errors. For the reconstructed 
covariance matrix 
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(10) 

the conservative lower bounds on homodyne detection 
efficiencies read 

rj A = 0.839, rj B = 0.780, rj c = 0.784. (11) 

On inserting the above explicit values into formula ([9| 
we obtain 
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(12) 

The separability of a Gaussian state with covariance ma- 
trix 7£ can be checked by calculating minimum symplec- 
tic eigenvalue jij of a covariance matrix corresponding 
to the density matrix of a partially transposed state p T i , 
j — A,B,C. For the matrix 7^ we obtain 

Ha = 0.85, hb = 1-07, ji c = 1-04. (13) 

We have jib > 1 and j±c > 1 which confirms the sepa- 
rability with respect to the splittings B\AC and C\AB. 



We should also check that the loss compensation does not 
lead to an unphysical covariance matrix. An explicit cal- 
culation confirms that 7^ satisfies the generalized Heisen- 
berg inequality 7^ + i J > and the lowest symplectic 
eigenvalue of 7^ reads jiq = 1-H > 1, hence jl is a 
physically allowed covariance matrix. 

If the true homodyne detection efficiencies rjj would 
be larger than the lower bounds rjj , then the true covari- 
ance matrix would be obtained by sending 7^ through 
local lossy channels with transmittances rjj /rjj. Since lo- 
cal losses cannot generate any entanglement, the B\AC 
and C\AB separability of the state is preserved for all 
actual homodyne detection efficiencies higher than the 
lower bounds rjj. 

Phase fluctuations 

Besides losses and inefficient detection, the homodyne 
detection can be also influenced by random phase fluctu- 
ations between signal and local-ocillator beams. A single- 
mode phase shift by angle <p couples the amplitude and 
phase quadratures X and P, 

X' = Xcos(/> + Psin(/>, P' = Pcos(/>-Xsin(/>. (14) 

We can conveniently describe this canonical transforma- 
tion by a 2 x 2 orthogonal matrix 

v J \ — sm0 cos0 J v ' 

multiplying the column quadrature vector (X, P) T . 

We shall consider independent random phase shift <pj 
for each mode, j = A, P, C. The three-mode phase shift 
operation is thus described by a block diagonal 6x6 
matrix 

5(0) = R^ A ) © R^b) © ij(0 c ), (16) 

where <fi = {(^a^c^d)- Let denote the covariance 
matrix in the absence of phase fluctuations. The phase 
shift S(4>) transforms the covariance matrix 7^ as fol- 
lows, 

i K = s{4>) lK s T {4>). (17) 

Since the mean values of all quadratures of the stud- 
ied state vanish, the observed covariance matrix 7^ can 
be determined by averaging the phase shifted covariance 
matrix 7^ over the random phase shifts 0, 

lL = {S{4>) lK S T {ct>))^. (18) 

As there is no reason for correlated phase fluctuations in 
our experiment, we shall assume that the random phase 
shifts are independent and mutually uncorrelated. To 
be specific, we shall assume that each phase shift obeys 
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Gaussian distribution with zero mean and variance a 2 . 
The relevant averages then read 
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(19) 
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Note that in the experimentally relevant limit of small 
phase fluctuations the exact shape of the probability dis- 
tribution becomes irrelevant and the above formulas may 
be considered as universally valid in the limit a <C 1. 



With the help of the expressions (19), (20), and (18) we 
obtain 



(i 



1L = e -V( 7K )+^— e - >-JV( lK )J T . 

(21) 

Here J denotes the three-mode symplectic form and 
^(jk) denotes a covariance matrix that is obtained from 
Ik by setting all covariances between quadratures of dif- 
ferent modes equal to zero. Hence V(jk) has a block 
diagonal structure with three 2x2 matrices on the main 
diagonal. The formula (21) can be inverted and after 
some algebra one finds 

o ( P <r 2 1 _ p 2a 2 

IK = e° 1L + 1 2 } P( 7 l) + — J— JV(hl) J T ■ 

(22) 

This expression allows us to undo the effect of phase fluc- 
tuations and determine the covariance matrix of a state 
before it was affected by the phase noise, c.f. Fig. S.l. 

If the state with covariance matrix 7^ is B\AC and 
C\AB separable, then the local phase fluctuations pre- 
serve this property because local single-mode operations 
cannot generate any entanglement from a separable state. 
In Fig. S.2 we plot the symplectic eigenvalues character- 
izing the entanglement properties of as a function of 
the strength of phase fluctuations a. We can see that a 
Gaussian state with covariance matrix 7^ remains B\AC 
and C\AB separable unless the phase fluctuations exceed 
the threshold value of <j t h ~ 7°. 

Note that the phase-diffused state with covariance ma- 
trix 7l can be slightly non-Gaussian due to phase fluc- 
tuations. Therefore, the covariance- matrix based separa- 
bility criteria valid for Gaussian states cannot be directly 
applied to this state (these criteria become only a nec- 
essary condition for separability but not a sufficient one 
for general non-Gaussian states). Instead, we have shown 
above that we can undo the effect of phase fluctuations 
on the reconstructed covariance matrix and determine 
the covariance matrix of a seed Gaussian state before the 
phase fluctuations. The separability properties of the 
phase diffused state can be then rigorously inferred by 
checking the separability of this seed Gaussian state. 



-2a z 




Figure 6: The lowest symplectic eigenvalue fio of covariance 
matrix 7^ (solid black line), symplectic eigenvalue \ia indi- 
cating entanglement in the A\BC splitting (black dot-dashed 
line) and the parameter min(/x J B,Mc) (black dashed line). We 
can see that the state remains B\AC and C\AB separable if 
the strength of random phase fluctuations is below a t h ~ 7° . 
The red dot-dashed line represents plot of min(/is, jic) for 
covariance matrix 70 with removed non-Gaussian fraction of 
noise added by random modulation in hot squeezing genera- 
tion. See text for details. 



Hot-squeezing generation 

Finally, let us consider the procedure for hot-squeezing 
generation by random displacements of the squeezed 
quadrature, see Fig. S.l. In the experiment, a great effort 
is made to achieve a Gaussian modulation and preserve 
the Gaussian nature of the state. Nevertheless, slight de- 
viations from perfectly Gaussian modulations may occur 
in practice and our goal here is to assess their impact. In 
the Heisenberg picture, the random modulation can be 
described by an addition of a noisy component xn to the 
squeezed amplitude quadrature of mode C, 



xc -> %c + x N . 



(23) 



The random variable x n has zero mean and let ajr denote 



its variance. We shall assume that the added noise xn 
contains both Gaussian and non-Gaussian part, 



x N 



Pg %n,g 



0" 



p G x N . 



(24) 



Here pc denotes the fraction of Gaussian modulated 
noise, xn,g is Gaussian random variable with zero mean 
variance and xn,g can exhibit an arbitrary distribu- 
tion with zero mean and variance aj^. 

The overall state generated by the noise addition 
becomes a statistical mixture of Gaussian states ob- 
tained by first adding the Gaussian part of the noise 
\/Pg%n,g and then randomly coherently displacing xc 
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by y/1 — pc %n- Note that the coherent displacements 
do not influence the entanglement properties of the state 
and can be set to zero by local unitary displacement op- 
erations. Therefore, if a three- mode Gaussian state ob- 
tained by considering only the Gaussian part of noisy 



Table I: Homodyne detector settings for complete reconstruc- 
tion of the three-mode covariance matrix. 



modulation in Eq. (24) is B\AC and C\AB separable, 
then also the total state including the full noisy modula- 
tion xn is B\AC and C\AB separable. 

Consider now the experimentally determined covari- 
ance matrix that already includes compensations for 
the inefficient homodyne detection and phase fluctua- 
tions. Since this matrix is by definition symmetric, it 
can be diagonalized by a suitable orthogonal transforma- 
tion O, 



7k = 0~/ d O 



(25) 



where 7^ is a diagonal matrix with eigenvalues of 7^ on 
the main diagonal. Given our method of preparation of 
the three-mode state by mixing uncorrelated single-mode 
states on balanced beam splitters, the six eigenvalues of 
7x should correspond to the variances of squeezed and 
anti-squeezed quadratures of the three input states, re- 
spectively. For instance, for a = 0.05 the six eigenvalues 
explicitly read 



Ai = 0.57, 
A 2 = 3.88, 



A 3 = 1.18, 
A 4 = 1.28, 



A 5 = 11.55, 
A 6 = 13.23. 



We can see that the eigenvalues have the expected struc- 
ture. Eigenvalues Ai and A2 correspond to the variances 
of squeezed and anti-squeezed quadratures of squeezed 
vacuum state in mode B. The variances A3 and A4 cor- 
respond to the vacuum state input in mode B. Ideally, 
we should have A3 = A4 = 1 for vacuum state but due 
to various experimental imperfections the eigenvalues are 
slightly larger. Finally, the last two eigenvalues A5 and A6 
correspond to the highly noisy hot-squeezed pseudother- 
mal state injected in mode C. 

We calculate the covariance matrix 70 of the effective 
Gaussian state obtained by considering only the Gaussian 
part xn,g of the random modulation xn by replacing Xj 
with pG^j where j = 5 or j = 6. Since it is difficult to 
reliably identify which of the two variances A5, A6 corre- 
sponds to the randomly modulated amplitude, we con- 
servatively rescale both these eigenvalues. Let 7^ denote 
the covariance matrix with A 5 and A6 replaced with pg^s 
andpcA6, respectively. Then 



7o = Oj D 7 



(26) 



We calculate the minimum symplectic eigenvalues 
and fic that characterize the separability properties of 
Gaussian state with covariance matrix 70 with respect 
to B\AC and C\AB splittings, respectively. The pa- 
rameter min(/i£, nc) is plotted as a red dot-dashed line 
in Fig. S.2. In this example we set pg = 0.75, i.e. 



setting # 


BHDa 


BHDs 


BHD C 


1 


X a 


X B 


X c 


2 


Pa 


Pb 


Pc 


3 


Pa 


X b 


X c 


4 


X a 


Pb 


X c 


5 


X a 


X b 


Pc 


6 


Xa+Pa 


X b +P b 


Xc+Pc 



we assume 25% of non-Gaussian modulation. We can 
see that our experiment is very robust against this non- 
Gaussian modulation and 25% of non-Gaussianity only 
very slightly reduces the maximum tolerable phase fluc- 
tuations. We estimate that in the actual experiment the 
no n- Gaussian fraction of the random modulation is much 
smaller than 25%. The dependences of the symplectic 
eigenvalues /io and Ha of 7o ° n the strength of phase noise 
(7 practically coincide with the results obtained for 
and therefore we do not plot them as separate curves in 
Fig. S.2. This coincidence is in full agreement with theo- 
retical predictions because both and /io are chiefly de- 
termined by the properties of the squeezed vacuum state 
in mode A and are not sensitive to change of fluctua- 
tions in input mode C provided that the noise in mode 
C remains high above the shot noise level. 



Statistical errors 

The three-mode covariance matrix jm is reconstructed 
from the quadrature measurements of the homodyne de- 
tectors. The quadrature operator measured by each ho- 
modyne detector Xj cos 0j + Pj sin 0j is determined by 
the relative phase 0j between local oscillator and signal 
beam. Amplitude quadrature Xj is measured for 6j = 
and phase quadrature Pj is recorded for 0j = tt/2. In 
order to probe the correlations between amplitude and 
phase quadratures of a given mode we also use the setting 
0j = 7r/4 which corresponds to measurement of Xj + Pj . 
In the experiment, we perform measurements for six dif- 
ferent settings of the three measured quadratures as in- 
dicated in Table S.I. These six settings are sufficient to 
determine all matrix elements of 7m- The variances and 
covariances of amplitude quadratures Xj are determined 
in the first step and this procedure is repeated also for 
the phase quadratures Pj in the second step. In the third 
(fourth, fifth) step we determine covariances of phase 
quadrature Pa (Pb, Pc) with the amplitude quadratures 
of the other two modes. Finally, the final sixth step yields 
the intramodal correlations (AXjAPj + APjAXj). A 
total number of 10 6 data was recorded for each measure- 
ment setting. 
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To determine the statistical errors of the symplectic 
eigenvalues fij calculated from the estimated covariance 
matrices, we have performed a Monte Carlo simulation 
of the whole experiment assuming Gaussian statistics of 
the measured quadratures. For each run of the Monte 
Carlo simulation we have reconstructed the covariance 
matrix and calculated the symplectic eigenvalues. This 
procedure was repeated 1000 times which provided a sta- 
tistical ensemble for each estimate jij. The mean values 
and statistical errors determined with the use of these 
ensembles read 

fi A = 0.849 ±0.001, 

fi B = 1.069 ±0.001, (27) 

lie = 1.036 ±0.001. 

This confirms that the statistical errors are very small 
compared to the deviations of jij from 1 hence the ob- 
served separability properties are statistically significant. 
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